#ifndef GEOMETRY_LINE_HPP
#define GEOMETRY_LINE_HPP

#include <cmath>
#include "point.hpp"

/* Represents a (directional) line that passes through two distinct points
 * in the 2-D plane, (x0, y0) and (x1, y1). The parametric equation is
 *     x = x0 + (x1 - x0) * t,
 *     y = y0 + (y1 - y0) * t.
 */
template <class T>
struct Line
{
    Point<T> P[2];
    Line() { }
    Line(const Point<T> &start_point, const Point<T> &end_point)
    {
        P[0] = start_point;
        P[1] = end_point;
    }

#if 0
    /// Returns the i-th end-point (i = 0, 1). 
    Point<T> & operator [] (size_t i) { return P[i]; }
    const Point<T> & operator [] (size_t i) const { return P[i]; }
#endif

#if 0
    /// ax + by + c = 0 
    T a() const { return P[0].y - P[1].y; }
    T b() const { return P[1].x - P[0].x; }
    T c() const { return P[0].x*(P[1].y-P[0].y) - P[0].y*(P[1].x-P[0].x; }
#endif

#if 0
    /* Checks whether a given parameter, t, is within bounds of the ray or
     * line segment. Returns one of the following values:
     *  0 : t is out of bounds
     *  1 : t is within bounds strictly
     *  2 : t is on the boundary
     */
    int hitTest(T t) const
    {
        switch (type)
        {
        case Ray:
            return (t > 0)? 1 : (t == 0)? 2 : 0;
        case Segment:
            return (t > 0 && t < 1)? 1 : (t == 0 || t == 1)? 2 : 0;
        default:
            return 1;
        }
    }
#endif

#if 0
    /// Checks whether the line is unbounded.
    bool isUnbounded() const 
    {
        return types[0]==EndPointUnbounded && types[1]==EndPointUnbounded;
    }

    /// Checks whether the line is a point or empty.
    bool isDegenerate() const 
    {
        return P[0] == P[1];
    }
#endif
};

/* Represents a (directional) line segment that passes through two distinct 
 * points in the 2-D plane, (x0, y0) and (x1, y1). The parametric equation is
 *     x = x0 + (x1 - x0) * t,
 *     y = y0 + (y1 - y0) * t.
 * where 0 <= t <= 1.
 */
template <class T>
struct LineSegment : public Line<T>
{
    LineSegment() { }
    LineSegment(const Point<T> &start_point, const Point<T> &end_point)
        : Line<T>(start_point, end_point) { }
};

#if 0
/// Represents a polygon on a plane.
template <class T, size_t Sides>
struct Polygon
{
    Point<T> P[Sides];
    size_t sides() const { return Sides; }
    Point<T>& operator[](size_t i) { return P[i]; }
    const Point<T>& operator[](size_t i) const { return P[i]; }
};
#endif

/* Returns the relative location of a point to a line, which is one of the
 * following values:
 *  1: if the point is on the left-hand side of the line, assuming the 
 *     direction of the line as the 'forward' direction.
 *  0: if the point is on the line or the line is degenerate.
 * -1: if the point is on the right-hand side of the line.
 */
template <class T>
int location(const Point<T> &P, const Line<T> &L)
{
    T a = L.P[0].y - L.P[1].y;
    T b = L.P[1].x - L.P[0].x;
    T c = -a * L.P[0].x - b * L.P[0].y;
    T det = a * P.x + b * P.y + c;
    return (det > 0)? 1 : (det < 0)? -1 : 0;
}

/* Returns the relative location of a point to a line segment, which is one 
 * of the following values:
 *  2: the point is on the extended line of the segment, but exceeds the
 *     segment's end-point (too far forward).
 *  1: the point is on the left-hand side of the line segment, assuming 
 *     the direction of the line as the 'forward' direction.
 *  0: the point is on the line segment or the line segment is degenerate.
 * -1: the point is on the right-hand side of the line segment.
 * -2: the point is on the extended line of the segment, but exceeds the
 *     segment's start-point (too far backward).
 */
template <class T>
int location(const Point<T> &P, const LineSegment<T> &L)
{
    int loc = location(P, (Line<T>)L);
    if (loc)
        return loc;
    T t = invMidPoint(P, L.P[0], L.P[1]);
    return (t > 1)? 2 : (t < 0)? -2 : 0;
}

#if 0
/* Gets the convexity of a quadrilateral. Returns one of the following:
 *  1: the quadrilateral is strictly vertex; that is, each pair of opposite
 *     vertices are on opposite sides of the line formed by connecting the 
 *     rest two vertices.
 *  0: the quadrilateral is degenerate; that is, at least three points are
 *     on the same line.
 * -1: the quadrilateral is strictly concave; that is, there exists a pair
 *     of opposite vertices that are on the same side of the line formed by
 *     connecting the rest two vertices.
 */
template <class T>
int convexity(const Polygon<T,4> &Q)
{
    Line<T> D0(Q.P[0], Q.P[2]);
    Line<T> D1(Q.P[1], Q.P[3]);
    
    int det = location(Q.P[1], D0) * location(Q.P[3], D0);
    if (det <= 0)
        return det;
    else
        return location(Q.P[0], D1) * location(Q.P[2], D1);
}
#endif

#if 0
/* Intersects two lines and returns the crossing point. If the lines do not
 * cross (i.e. they are parallel), returns 0. If the lines are identical,
 * returns 2 and no crossing point is stored.
 */
template <class T>
int intersect(const Line<T> &L1, const Line<T> &L2, Point P[1])
{
    T a1 = L1.a, b1 = L1.b, c1 = L1.c;
    T a2 = L2.a, b2 = L2.b, c2 = L2.c;
    T det = a1*b2 - a2*b1;
    T p = b2*c1 - b1*c2;
    T q = a2*c1 - a1*c2;

    if (det == 0) /* identical or parallel */
        return (p == 0 && q == 0)? 2 : 0;

    P[0].x = p / det;
    P[0].y = q / det;
    return 1;
}
#endif

#if 0
/* Intersects two lines and stores the unique crossing point if one is found.
 * Otherwise, no crossing point is stored. Returns one of the following:
 *  0: no crossing point is found, either because the lines are parallel, or 
 *     because the crossing point would exceed the boundary of one of the
 *     line segments. Also, if any of the lines is degenerate, returns 0.
 *  1: a unique crossing point is found and stored in P.
 *  2: there are infinitely many crossing points because the lines overlap.
 */
template <class T>
int intersect(const Line<T> &L1, const Line<T> &L2, Point<T> P[1])
{
    /* Return failure if L1 or L2 is degenerate. */
    if (L1.P[0] == L1.P[1] || L2.P[0] == L2.P[1])
        return 0;

    T x1 = L1[0].x, y1 = L1[0].y, p1 = L1[1].x-L1[0].x, q1=L1[1].y-L1[0].y;
    T x2 = L2[0].x, y2 = L2[0].y, p2 = L2[1].x-L2[0].x, q2=L2[1].y-L2[0].y;
    T det = p2*q1 - p1*q2;
    T d1 = (x1-x2)*q2 - (y2-y1)*p2;
    T d2 = (y2-y1)*p1 - (x2-x1)*q1;

    if (det != 0) /* non-parallel or overlapping */
    {
        T u = d1 / det;
        T v = d2 / det;
        if (L1.hitTest(u) && L2.hitTest(v))
        {
            P = midPoint(L1[0], L1[1], u);
            return 1;
        }
        return 0;
    }

    if (d1 != 0 || d2 != 0) /* no solution because the lines are parallel */
        return 0;

    /* Now that L1 and L2 are overlapping, check the bounds. */
    if (L1.type == LT_DEFAULT || L2.type == LT_DEFAULT)
        return 2;

    /* Find the parameter of the end points of line 2 on line 1. */
    T u1 = invMidPoint(L2.P[0], L1.P[0], L1.P[1]);
    T u2 = invMidPoint(L2.P[1], L1.P[0], L1.P[1]);
    if (u2 >= u1)   /* L1 and L2 are same direction */
    {
        if (L1.type == LT_RAY && L2.type == LT_RAY)
        {
        }
    }
    else            /* L1 and L2 are opposite direction */
    {
    }

    if (L1.type == LT_SEGMENT && L2.type == LT_SEGMENT)
    {

    }

}
#endif

/* Intersects two line segments and stores the unique crossing point if one
 * is found. Returns one of the following values:
 *  0: no crossing point is found, either because the lines are parallel, or 
 *     because the crossing point would exceed the boundary of at least one 
 *     of the line segments.
 *  1: a unique crossing point is found and stored in P.
 *  2: there are infinitely many crossing points because the line segments
 *     (partially) overlap.
 */
template <class T>
int intersect(const LineSegment<T> &L1, const LineSegment<T> &L2, Point<T> P[1])
{
    /* Handle the case where one of the line segments is degenerate. */
    if (L1.P[0] == L1.P[1])
        return location(P[0] = L1.P[0], L2) == 0 ? 1 : 0;
    if (L2.P[0] == L2.P[1])
        return location(P[0] = L2.P[0], L1) == 0 ? 1 : 0;

    T x1 = L1.P[0].x, y1 = L1.P[0].y, p1 = L1.P[1].x-L1.P[0].x, q1=L1.P[1].y-L1.P[0].y;
    T x2 = L2.P[0].x, y2 = L2.P[0].y, p2 = L2.P[1].x-L2.P[0].x, q2=L2.P[1].y-L2.P[0].y;
    T det = p2*q1 - p1*q2;
    T d1 = (x1-x2)*q2 - (y2-y1)*p2;
    T d2 = (y2-y1)*p1 - (x2-x1)*q1;

    if (det != 0) /* non-parallel or overlapping */
    {
        T u = d1 / det;
        T v = d2 / det;
        if (u >= 0 && u <= 1 && v >= 0 && v <= 1)
        {
            P[0] = midPoint(L1.P[0], L1.P[1], u);
            return 1;
        }
        return 0;
    }

    if (d1 != 0 || d2 != 0) /* no solution because the lines are parallel */
        return 0;

    /* Find the parameter of the end points of line 2 on line 1. */
    T u1 = invMidPoint(L2.P[0], L1.P[0], L1.P[1]);
    T u2 = invMidPoint(L2.P[1], L1.P[0], L1.P[1]);
    if (u1 <= u2)
    {
        if (u1 < 1 || u2 > 0)   /* infinite */
            return 2;
        if (u1 > 1 || u2 < 0)   /* none */
            return 0;
        if (u1 == 1)
            P[0] = L2.P[0];
        else
            P[0] = L2.P[1];
        return 1;               /* unique */
    }
    else                        /* u1 > u2 */
    {
        if (u2 < 1 || u1 > 0)   /* infinite */
            return 2;
        if (u2 > 1 || u1 < 0)   /* none */
            return 0;
        if (u2 == 1)
            P[0] = L2.P[1];
        else
            P[0] = L2.P[0];
        return 1;               /* unique */
    }
}

#if 0
/* Finds the locus of a point P such that |PL1|/|PL2| = p/q. */
int splitAngle(Line result[2], const Line &L1, const Line &L2, double p, double q)
{
    double a1 = L1.a, b1 = L1.b, c1 = L1.c;
    double a2 = L2.a, b2 = L2.b, c2 = L2.c;
    double u = p * std::sqrt(a1*a1+b1*b1);
    double v = q * std::sqrt(a2*a2+b2*b2);

    result[0] = Line(v*a1-u*a2, v*b1-u*b2, v*c1-u*c2);
    result[1] = Line(v*a1+u*a2, v*b1+u*b2, v*c1+u*c2);
    return 2;
}

/* Solves the linear system of equations:
 *   a11*x + a12*y = b1
 *   a21*x + a22*y = b2
 * Returns the determinant of the system.
 */
double solveLinearSystem(
    double *x, double *y,
    double a11, double a12, double b1,
    double a21, double a22, double b2)
{
    double det = a11*a22 - a12*a21;
    *x = (a22*b1 - a12*b2) / det;
    *y = (a11*b2 - a21*b1) / det;
    return det;
}

int intersect(Point P[1], const Line &L1, const Line &L2)
{
    solveLinearSystem(&P[0].x, &P[0].y, L1.a, L1.b, -L1.c, L2.a, L2.b, -L2.c);
    return 1;
}

/// Gets a line that is perpendicular to P1P2 and passes through P0.
Line getPerpendicular(const Point &P1, const Point &P2, const Point &P0)
{
    // The equation is: (x2-x1)(x-x0) + (y-y0)(y2-y1) = 0.
    Line L(P2.x-P1.x, P2.y-P1.y, (P1.x-P2.x)*P0.x + (P1.y-P2.y)*P0.y);
    return L;
}
#endif

#if 0
Line transform(const Line &L, const Transformation &t, const Transformation &tinv)
{
    double A = L.a, B = L.b, C = L.c;
    double a11=tinv.a[0][0], a12=tinv.a[0][1], a13=tinv.a[0][2];
    double a21=tinv.a[1][0], a22=tinv.a[1][1], a23=tinv.a[1][2];
    Line LL(A*a11+B*a21, A*a12+B*a22, A*a13+B*a23+C);
    return LL;
}
#endif

#endif /* GEOMETRY_LINE_HPP */
